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The  Generation  of  Three-Dimensional  Body-Fitted 
Coordinate  Systems  For  Viscous  Flow  Problems 

By 

Z.  U.  A.  Wars!  ,  C.  W.  Mastin  and  J.  F.  Thompson 


Abstract 

The  proposed  set  of  equations  (refer  to  the  enclosed  papers  for 
detail)  which  generate  a  series  of  surfaces  between  a  given  inner  and 
outer  body,  have  been  programed  on  CRAY-1.  An  extensive  program  test¬ 
ing  has  been  carried  out  to  make  the  program  usable  for  general  body 
shapes.  For  example,  two  methods  have  been  developed  to  establish  the 
correspondence  between  the  points  of  the  inner  and  outer  surfaces,  a 
method  has  been  developed  to  find  the  first  partial  derivative  of 
x,  y,  z  with  respect  to  the  coordinate  along  the  surface.  All  these 
methods  are  based  on  sound  mathematical  basis  and  have  not  been  chosen 
arbitrarily.  Work  has  been  started  on  complicated  multibody  config¬ 
urations,  such  as  the  wing-body  combination.  Here  the  interfacing  of 
coordinates  having  sufficient  derivative  smoothness  is  the  most  impor¬ 
tant  problem. 

In  the  period  of  this  report,  a  thorough  analysis  on  all  sorts  of 
mathematical  models  for  coordinate  generation  has  been  completed.  This 
analysis  uncovers  those  differential  relations  which  must  invariably 
be  satisfied  by  the  metric  coefficients  no  matter  which  method  is  used  - - J 
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to  generate  the  coordinates.  The  final  computational  code  to  be 
developed  in  the  current  year  of  effort  is  expected  to  reflect  all 
these  achievements  for  practical  utilization.  ~ 


Previously  reported  examples  of  surface  grid  construction  have 
dealt  with  simple  quadratic  surfaces.  To  test  the  versitility  of  the 
method,  our  program  has  been  modified  to  accept  the  surface  data 
generated  by  Craidon  [1]^  Aside  from  the  more  complex  computational 
region,  which  gave  rise  to  additional  coordinate  singularities,  the 
grid  generation  procedure  was  unchanged.  That  is,  the  grid  was 
generated  using  cubic  splines  with  an  elliptic  system  used  for  smoothing. 
Several  coordinate  surfaces  for  the  region  about  a  spline  generated 
wing-fuselage  configuration  have  been  plotted.  A  plot  of  that  con¬ 
figuration  is  given  in  Figure  1(a).  Figure  1(b)  illustrates  the 
continuation  of  coordinate  lines  from  the  body  to  the  outer  boundary. 

Views  of  coordinate  surfaces  from  the  upstream  direction  are  also  included 
Figure  1(c)  is  a  surface  surrounding  the  aircraft,  and  Figure  1(d)  is 
a  surface  intersecting  the  trailing  edge  of  the  wing. 


+ 
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51.  INTRODUCTION 

This  paper  examines  in  detail  the  analytical  aspects  of  three  distinct 
methods  of  coordinate  generation  based  on  partial  differential  equations, 
in  either  two  or  three  dimensions.  The  first  method  is  based  on  the  Gauss 
equations  of  a  surface  under  the  constraint  of  the  Beltrami's  second  order 
equations.  These  equations  have  been  structured  in  such  a  way  that  an 
automatic  connection  is  established  between  the  succeeding  generated  surfaces. 
The  second  method  is  a  re-examination  of  those  equations  which  are  based 
on  the  inhomogeneous  Laplace  equations.  This  analysis  reveals  a  new  form 
for  the  terms  which  play  a  role  in  the  concentration  of  coordinate  lines 
and  in  the  adaptive  coordinate  system  generation.  The  third  method  pertains 
to  a  set  of  equations  in  the  metric  coefficients  which  is  obtained  by  setting 
the  Riemann* s  curvature  tensor  to  zero. 

The  problem  of  generating  spatial  coordinates  by  numerical  methods  is  a 
problem  of  much  interest  in  practically  all  branches  of  engineering  and 
physics.  At  present  a  number  of  techniques  are  under  active  development  for 
the  generation  of  two  and  three-dimensional  coordinates  in  the  regions 
between  two  or  a  number  of  arbitrary  shaped  bodies.  Among  these  efforts 
two  easily  discernable  groups  can  be  formed,  (i)  the  methods  based  on 
elliptic  PDE's,  and  (ii)  algebraic  methods.  In  the  first  group,  a  set  of 
inhomogeneous  Laplace  equations  is  taken  as  the  basic  generating  system. 

These  equations  are  then  inverted  and  solved  for  the  Cartesian  coordinates. 

^Professor 


2 


, ,  Some  very  useful  results  based  on  this  line  of  approach  started  with  the 

•  1  2 
work  of  Winslow  ,  have  been  obtained  by  Thompson,  et  al.  (TTM  method), 

Steger,  et  al.3,  Yu4,  Graves^,  and  Thomas  .  For  an  extensive  bibliography 

7 

refer  to  Thompson,  et  al.  In  the  second  group  of  methods,  the  grid  points 
in  space  are  generated  by  interpolating  and  blending  functions  starting 

from  the  given  boundary  data.  This  line  of  approach  has  been  followed  by 

8  9  10 

Eiseman  ,  Smith,  et  al.  ,  Erickson  ,  and  others. 

In  this  paper  we  consider  only  the  analytical  aspects  of  the  differential 

equation's  approach  to  coordinate  generation.  The  main  effort  here  is 

to  present  only  those  results  which  are  of  permanent  interest  to  the  workers 

in  the  field  of  coordinate  generation.  The  proposed  equations  in  any  one 

of  the  groups  have  not  been  arbitrarily  selected  to  generate  some  sort  of 

coordinates.  These  equations  are  in  fact  those  which  every  numerically 

or  analytically  generated  coordinates  must  satisfy.  The  reader  will 

find  that  some  large  portions  of  sections  3  and  5  have  new  results  and 

11  12 

are  based  on  the  work  by  Warsi  '  .  In  sections  3  and  5  a  number  of 

exact  solutions  have  been  obtained  which  can  be  used  to  provide  a  testing 
ground  for  different  numerical  schemes. 


§2.  NOTATION  AND  BASIC  FORMULAS 

In  this  paper  any  general  curvilinear  coordinate  system  will  be  denoted 
by  a  superscript  index  notation,  such  as  x1.  However,  when  an  expression 
has  been  expanded  out  in  full  and  there  is  no  need  for  an  index  notation  then 
we  shall  use  the  symbols 


1 


x 


=  £ 


3 

n  ,  x  =  5  • 


The  rectangular  Cartesian  coordinates  (x,  y,  z)  which  determine  the  position 
vector  r,  i.e., 

r  =  r(x,y,z) 


will  be  denoted  by  the  subscripted  variable  x.,  where  x  =  x,  x  »  y, 

1  1  ^ 

x3  *  z. 

Two  similar  indices,  one  appearing  as  a  subscript  and  the  other  as  a 

superscript  will  always  imply  summation  over  the  range  of  index  values;  e.g., 

A,  .B1  «  A,  .B1  +  A..B2  +  A,.B3. 
ij  lj  2j  3 j  ^  3 

In  an  Euclidean  space  (E  or  E  ) ,  the  covariant  base  vectors  a^  are  given 

by 


so  that 


*1  “  '  *2  =  '  -3  ~  fc  ' 


(lb) 


where  a  variable  subscript  will  denote  a  partial  derivative.  Using  the 
Riemannian  metric,  the  formula  for  the  length  element  ds  is  given  by 

(ds)2  =  gi^dx1dx^  , 

where,  because  of  the  Euclidean  nature  of  the  space  the  metric  coefficients 
are  given  by 


3r  9r 

g.  .  =  a .  *a  .  =  ~~r  ’  — . 

13  -1  '3  3X1  ax3 


(2a) 


The  coefficients  q  =  g..  are  the  covariant  components  of  the  metric  tensor. 
13  31  ii 

The  contravariant  components  g  are  related  with  g  .  through  the  equation 

iD 


=  ,3 

9  gik  6k 


(2b) 


where  6  3  (the  Kronecker  deltas)  are  the  mixed  components  of  the  metric  tensor. 
Using  (2b)  we  define  the  contravariant  base  vectors  as 


l  13 
a  =  g  a . 

O 


(2c) 


The  quantities  g  and  g  defined  as 


g  «  det(g_)  , 


g  *  det (g13 )  , 


(3a) 

(3b) 


are  related  as 


(4) 


For  a  three-dimensional  space 

9  "  911922933  +  2g12913923  "  (g23)2gll  ‘  (g13)2g22  "  <gl2)2g33*  (5) 


4 


Introducing  the  quantities, 


G1  =  922933 


(923)  ' 


G2  =  911933 


(913>  ' 


G3  911922 


(912}  ' 


G4  =  913923  "  912933  ' 


G5  912923  ”  913922  ' 


G6  912913  ~  923911  ' 


we  have,  on  solving  Eqs.  (2b) , 


11  „  .  22  33 

g  =  G  /g  ,  g  =  G2/g  ,  g  =  G3/g  , 


g12  =  G./g  ,  g13  =  G  /q  ,  g23  *  G  /g 


6/y  ’ 


The  space  Christoffel  symbols  of  the  first  and  second  kind  respectively 


are  given  by 


Using  (8b) ,  we  have 


1  3gik  39ik  3gii 

[ij,k]  =  j(-^f  +  -f1  -  -f1)  , 

2  3x3  9xX  3x* 


k£.r  .  .  .  , 

Fij  =  9 


3a ,  . 

.  r* ,  a  . 


In  the  case  of  a  two-dimensional  surface  embedded  in  a  three-dimensional 
space,  we  shall  use  the  Greek  indices  a,  8,  etc.  (with  the  exception  of  v) 
with  the  stipulation  that  they  assume  only  two  values.  Thus  the  surface 
Christoffel  symbols  of  the  first  and  second  kind  are  respectively  given  by 

.  ..  I.9ga5  3g66  3gaB.  ,*n> 

la8,6j  *  r-( — tt  +  r - 7>  •  iaa' 

2  3x°  3x°  3x6 


5 


Tag  “  960ta6  r  Ob) 

where  for  the  purpose  of  clarification  we  have  used  the  symbol  T  (upsilon) 
to  denote  the  surface  Christoffel  symbols  of  the  second  kind  in  (9b)  and 
net  by  T  as  in  (8b) . 

In  the  process  of  formulation  of  a  3D  coordinate  generation  problem,  it 

is  helpful  to  imagine  the  coordinates  of  a  point  in  space  as  the  intersection 

of  three  distinct  surfaces  on  each  of  which  one  coordinate  is  held  fixed. 

12  3 

Using  the  convention  of  a  right-handed  coordinate  system  x  ,  x  ,  x  or 

£,  n»  5,  we  introduce  the  notation  2 ^  as  a  surface  on  which  the  coordinate 
v 

x  =  const.,  such  that 
• 

2  3 

v  =  1  implies  that  (x  ,x  )  are  in  the  surface, 

v  =  2  implies  that  (x3,  x^)  are  in  the  surface, 

v  =  3  implies  that  (x1,  x2)  are  in  the  surface. 

Thus,  the  unit  normal  vector  on  the  surface  2 ^  is 

n(v)  =  (ra*rg)/|raxre|  ,  (10) 

where 

v  =  l:a=2,{5=3  (surface  x1  =  £  =  const.)  , 

2 

v  =  2  :  ct  =  3  ,  8=1  (surface  x  =  n  =  const.)  ,  (11) 

v=3:a=l,8=2  (surface  x3  =  £  =  const.)  . 

All  other  quantities  and  formulas  which  appear  in  the  rest  of  the  paper  have 

12  13 

been  defined  where  they  first  appear.  Refer  also  to  Harsi  and  Eisenhart 

53.  GENERATING  DIFFERENTIAL  EQUATIONS  BASED  ON  GAUSS  EQUATIONS 

In  this  section  our  aim  is  to  develop  a  method**  for  the  generation  of 
3D  coordinates  wherein  a  series  of  surfaces  are  generated  on  each  of  which 
two  previously  designated  coordinates  vary  while  the  third  coordinate  remains 
fixed.  This  method  must  also  be  structured  in  such  a  way  that  the  variation 


of  the  third  coordinate  from  one  generated  surface  to  the  next  is  fully 
reflected  in  the  system  of  generating  equations.  With  this  aim,  we  start 
from  the  equations  of  Gauss^'^  which  for  a  surface  xV  =  const.,  are  given  by 


r  =  r.  +  b  n^ 
.ag  o0  .6  a0. 


where  the  variations  of  o,  0  and  the  range  of  6  with  v  follows  the  scheme  in 

(11) .  The  quantities  b  are  the  coefficients  of  the  second  fundamental 

“6  v  (v) 

form  of  the  surface .  Since  on  the  surface  x  =  const . ,  the  vector  n  is  • 


orthogonal  to  the  surface  vectors  r^  ,  hence 


b  =  n  ^  *r  .  (13) 

a0  .  .a0 

To  fix  ideas,  we  envisage  a  surface  which  is  formed  of  the  coordinate  lines 
4,  n  and  on  which  t,  =  const.  Dropping  the  index  v,  Eq.  (12)  yields  the 
three  equations 


=  tii^6  +  s"  ' 


(14a) 


Un  =  T!2-<5  +  70 


(14b) 


r  =  T  ,r.  +  Un 
.nn  22.6 


where  the  index  6  now  varies  from  1  to  2,  and 


(14c) 


S  =  bU  ,  T  =  bi2  *  u  =  b22  ’ 


Here  n  is  orthogonal  to  both  r^  and  r^,  and  the  coefficients  of  the  first 
fundamental  form  of  the  surface  are  g^ ,  g^2»  and  g22;  each  evaluated  at 


C  *  const.  Obviously 


<>11  -  '  912  ■  “{VWVn  '  *22  ■  vv\  • 


If  Eqs.  (14)  are  considered  as  the  first  order  partial  differential 
equations  in  rr  and  r  ,  then  we  must  also  consider  the  Weingarten  equations 


If  now  g^,  g12,  g22,  bn,  b12,  b22  are  arbitrarily  prescribed  then  the  set  of 
Eqs.  (14)  and  (17),  which  represent  fifteen  scalar  equations  for  the  nine 
scalars  (r^  ,  r^  ,  n ) ,  form  an  overdetermined  system.  Consequently  one 
has  to  impose  the  compatability  requirements 

(r  )  =  (r  ) 

-a6  y  ~aY  0 

for  all  values  of  a,  £,  y  from  1  to  2.  This  operation  leads  to  the  Mainardi- 
Codazzi  equations  ^nd  the  theorems  egregium  of  Gauss  which  are  higher  order 
equations  and  are  not  very  suitable  for  the  purpose  of  numerical  solution . 

We  therefore  return  to  the  Gauss  equations  (14)  and  ask  the  question:  Is 
it  possible  to  develop  a  method  which  centers  around  the  Gauss  equations  and 
is  simple  to  implement  numerically?  The  answer  is  in  affirmative  if  we 
manipulate  Eqs.  (14)  as  follows. 

Multiplying  Eq.  (14a)  by  g22'  E9*  by  ~2^i2  and  Eq‘  by  9j_i 

and  adding  the  three  equations,  we  get 

£r  -  -(<42{>r{  ♦  <42n>*;„l  G3 

+  <g22S  -  2gl2T  .  g11u)n  ,  (16) 

where  £,  is  the  second  order  differential  operator, 

g223u  "  2g1235n  +  gH9nn  ' 


and  A2  is  the  second  order  differential  operator  of  Beltrami. 
xV  =  const.,  (refer  to  the  scheme  in  (ID), 


For  any  surface 


*2  l3a{/5"(9663°  " 


+  M—  (g  90  -  g  ))] 

t  aa  B  aB  ot 

v 


(19a) 


In  particular  for  the  surface  <;  *  const,  we  drop  the  enclosed  superscript  and 


8 


where 


9 


\  =  n*r^  =  Xx^  +  Yy^  +  Zz^  , 

*  ‘  (Ytzn  '  Vt1^  • 

Y  =  (x  z  -  x,z  )//g7  , 
n  £  £  n  3 

z  -  <*{yn  -  V?)/^  • 


(23) 


(24) 


Thus,  by  using  the  forms  in  (22)  we  have  established  a  connection  with  the 
coordinate  £  which  changes  from  one  surface  to  the  next.  We  now  rewrite 
Eq.  (18)  as 


and  take  them  as  the  fundamental  generating  equations  for  the  coordinates  in 

a  surface.  It  must  be  noted  that  is  not  a  2D  Laplace  operator  except 

when  the  surface  degenerates  into  a  plane  having  no  dependence  on  z. 

It  is  a  well  known  result  in  differential  geometry  that  the  isothermic 

coordinates  in  a  surface  satisfy  Eqs.  (27)  identically.  The  isothermic 

coordinates  £  and  n  are  those  orthogonal  coordinates  in  a  surface  which  yield 

g  _  *  g  .  The  situation  here  is  parallel  to  the  choice  of  the  Laplace 
22  2  2 

equations  V  0,  7  n  *  0  for  the  generation  of  plane  curvilinear  coordinates. 


(e.g. ,  the  TTM  method  ),  which  are  also  satisfied  identically  by  the 
conformal  coordinates  in  a  plane.  This  does  not  mean  that  the  Laplace 
equations  are  suitable  only  for  the  generation  of  conformal  coordinates. 

In  fact,  as  is  evidenced  by  the  available  body  of  numerical  results,  the 
Laplace  equations  are  capable  of  generating  very  general  coordinates  in 
arbitrary  domains.  Therefore,  there  looks  to  be  no  apparent  reason  why 
Eqs.  (27)  should  not  form  the  basic  generating  system  for  general  coordinates 
in  a  surface.  The  analytical  solutions  given  in  this  paper  and  the  numerical 
results  given  in  Warsi  and  Ziebarth"^  support  this  contention. 

Having  chosen  Eqs.  (27)  as  the  generating  system,  the  equation  for  the 
determination  of  the  Cartesian  coordinates,  viz.,  Eq.  (25),  becomes 

£r  =  nR  .  (28) 

The  three  scalar  equations  in  expanded  form  are 


22  a 


(29a) 


g22^C  ‘  2gi2y£o  +  gllynn  =  YR 


(29b) 


g22zU  '  2g12z£n  +  gllznn  =  ZR  ' 


(29c) 


where  X,  Y,  Z,  and  R  have  been  defined  in  Eqs.  (24)  and  (26).  It  must  be 
noted  that  by  cyclic  permutations,  equations  similar  to  Eqs.  (29)  can  be 
written  for  the  surfaces  rj  =  const,  and  £  =  const.  However,  only  one  set, 
e.g.,  Eqs.  (29),  is  sufficient  provided  that  we  are  able  to  take  care  of  the 
derivatives  r^  appearing  in  R. 

The  set  of  Eqs.  (29)  form  a  consistent  set  of  equations  for  the  deter¬ 
mination  of  x,  y,  z  under  the  prescribed  boundary  conditions.*  For  an 

analytical  understanding  of  these  equations  we  open  the  differentiations  of 

3  3  3 

the  metric  coefficients  in  the  formulae  for  r^,  T and  ?22.  Thus 


ril  =  + 


6y«  +  Y2« 


(30a) 


+  By 


4  YZ 


Cn  ' 


(30b) 


•Refer  to  conment  (i)  at  the  end  of  the  paper. 
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'22 


coc 


nn 


+  6 ynn  +  yz 
nn 


nn 


(30C) 


where 


a  =  (G5x^  +  G6xn  +  G^x^J/g  , 

6  =  (G5y£  +  Vn  +  G3V/g  ' 


Y  =  (G^Zfr  +  G6zn  4  G3zc)/g  . 


Substituting  Eqs.  (30)  in  (26)  and  after  arranging  the  terms  we  can  rewrite 
Eqs.  (29)  as  a  quasilinear  system, 


..  32x. 

&  _ 1-  = 


a$  3xa3x6 


=  0  ,  i  =  1,  2,  3 


(31) 


where  x  =  x,  x  =  y,  x  =  z  and  there  is  an  implicit  sum  on  j  from  1  to  3 

123  j* 

and  on  a,  6  from  1  to  2.  The  coefficients  A^g  depend  on  the  metric  coefficients 
gll'  g12'  g22  and  °n  those  geometric  quantities  which  depend  only  on  the  first 
partial  derivatives.  For  example 

A11  *  g22(l-alX)  ,  A11  *=  -2g12(l-aAX)  ,  etc.,  etc. 


Equations  (31)  are  three  equations  in  three  unknowns  with  two  independent 
variables.  Refer  to  Petrovsky16  for  the  classifications  of  such  equations. 


S3. 2  Coordinate  redistribution  (concentration) 

Before  discussing  the  basic  solution  algorithm  for  the  set  of  Eqs.  (29) 
it  is  important  to  study  the  effect  of  a  coordinate  transformation  which 
produces  a  nonuniform  distribution  of  coordinates.  Again  using  indexed 
quantities,  let  xa  be  another  coordinate  system  defined  as 

-a  _a .  12 

x  =  x(x,x),a  =  l,2, 

with 

axa 

det(-22~)  *  0  . 

9x6 

Using  to  mean  either  x,  y,  or  z,  we  have 


3x '  3x 


a2x. 


3  Xi  35s  3xY  axi  32xY 
+ 


3xa3x^  3x^3xY  3x°  3x^  3xY  3xQ3x 


Also, 


a6  ,00  3x°  3x^ 

g  -  9  — a  — -  • 

3x  3x 


Now,  Eqs.  (29)  can  be  written  in  a  compact  form  as 


,2 

AS  8  ’i  .  *  „ 

9  I - a  =  ~  X,  * 


3x°3xe  S  1 


where 


X1  =  X  ,  x2  *  Y  ,  X3  =  Z 


and 


11  12  22 
9  =  g22/G3  ,  g  =  -g12/G3  ,  g 


On  coordinate  transformation  we  have 


G3  -  G3/(I?)  ,  R  -  R/(D)Z  ,  Xi 


where 


J-  3x^  3xi  3x^  3x^_ 

,-l  _2  "  _2  l  • 

3x  3x  3x  3x 


Thus  Eq.  (33)  becomes 


a  ^xi  3*.  c 

g“B  - |  ♦  5"°  PY  -i  -  -S 

3Jt°3xB  3xY  g3 


y  3x“  3xb  a2«T 


te"  3X0  3x“3x6 


where 


Using  equations  similar  to  (34)  in  the  new  coordinate  system,  Eq.  (36)  yields 
the  equations 


£x  =  XR  , 

(38a) 

£y  -  yr  , 

(38b) 

£z  =  ZR  , 

(38c) 

where 


£ 


g22  hi  "  2gi23£n  + 


gll3nn 


+  P3-  +  Q9- 

£  n 


(39) 


P  = 


922P11 


"  2g12P12  +  gllP22 


(40a) 


5  *  ^22^11  '  2512P?2  *  5up22  •  (40b> 

and  X,  Y,  2,  and  R  have  exactly  the  same  expressions  as  in  (24)  and  (26a) 
in  the  new  coordinate  system. 

y 

The  structure  of  the  terms  is  quite  revealing  particularly  in  those 

situations  when  it  is  desired  to  redistribute  an  already  existing  coordinate 
a 

system  x  so  as  to  achieve  a  desired  concentration  or  expansion  of  the 
coordinates  x  .  Though  still  a  forcing  function  behavior  for  P^o  has  to 
be  prescribed,  the  user  is  at  least  aware  of  its  structure,  that  is,  it  must 
be  composed  of  the  product  of  two  first  partial  derivatives  and  a  second 
partial  derivative.  These  considerations  may  be  important  in  the  adaptive 
coordinate  systems.  In  other  cases  P^  may  be  prescribed  arbitrarily.  One 
such  case  has  been  treated  numerically  in  Ref.  15.  (Refer  also  to  S3.) 


S3. 3  Morphology  of  A  Solution  Algorithm 

The  discussion  that  follows  pertains  to  the  case  when  it  is  desired  to 
generate  the  3D  curvilinear  coordinates  between  two  artibrary  shaped  smooth 
surfaces.  As  is  shown  in  Fig.  1,  let  the  surface  coordinates  of  the  inner 
body  n  ■  n_.  and  of  the  outer  body  n  *  n  be  the  same  coordinates.  Because 

O  00  1 

of  the  right-handedness  of  the  coordinate  triple  (£,n,C) ,  the  ordered  pair 
(<;,£)  is  taken  as  a  positive  ordered  pair  on  both  the  surfaces.  Since  both 
the  surfaces  n  *  nn  and  n  *  n  are  known  either  analytically  or  numerically, 
so  that 


n 


« 


R 


'B 
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r  *  r  (4,0  ;  n  -  n  :  r  »  r  (4,4)  , 

•.  00  ^  ^®D 


(41) 


and  hence  the  needed  partial  derivatives  with  respect  to  4  and  4  are  directly 
available  at  the  surface^ 


Figure  1.  Selection  of  coordinates  on  the  inner  and  outer  boundaries. 

For  the  computation  of  r  in  the  field  one  must  first  note  that  the 

(2) 

coordinate  4  may  not,  in  general,  satisfy  the  Beltrami's  equation  A  2  4 
Consequently,  r^  must  satisfy  the  equation 


=  0. 


(2)r  +  G2(Af}4)r5 


G2(k1(2)  +  k22))n(2) 


From  this  equation  we  devise  a  weighted  integral  formula 


u 


lfl(n)(!44)B  +  f2(n)<?44)-ld« 


(42a) 


where 


[-^  <*«>  *x'2>>n‘2> 
gu  1  2  -  gu  -« 


^33 

*11  -« 


/ST 


a  *11, 


3  ,913, 


and 


gn 


(42b) 


fl(V 


1  , 


0  •  VV 


(42c) 

Referring  to  Fig.  2(a),  we  now  solve  Eqs.  (29)  or  (38)  for  each  4  *  const.. 


0  ,  f, (n_)  *  1 
2  00 


by  prescribing  the  values  of  x ,  y  and  z  on  the  lower  curve  and  the  upper 

curve  which  represent  the  curves  on  B  and  •  respectively.  In  Fig.  2(b) 

and  are  the  cut  lines  on  which  periodic  boundary  conditions  are  to  be  imposed. 


r 

*> 


1 


-J 


Figure  2.(a)  Topology  of  the  given  surfaces,  (b)  Surface  to  be  generated. 


§3.4  Exact  solutions 

The  following  two  examples  demonstrate  that  the  proposed  set  of  generating 
equations  (27)  or  equivalently  the  set  of  equations  (29)  or  (38)  are  consistent 
and  provide  nontrivial  solutions. 

Example  1:  Isothermic  coordinates  on  a  unit  sphere. 

Let  the  surface  coordinates  of  a  unit  sphere  be  denoted  as  £,  £,  where 
the  order  (?;,£)  forms  a  right-handed  system.  Since  our  objective  is  to 
provide  isothermic  coordinates  which  are  orthogonal,  we  assume 

X  »  <MC)  ,  y  =  f(C)cos  C  ,  Z  ■  f  (C) sin  £  ,  (43a) 

so  that 

f2  +  1?  -  1  .  (43b) 

Calculating  the  metric  coefficients  and  the  surface  Christoffel  symbols 
based  on  the  assumed  form  (43a) ,  we  find  that  the  equations  *  0  and 
£> *  0  are  satisfied  provided  that 
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Eliminating  <J*  between  (43b, c)  ,  we  get 


2  2  2 
f*  =  (1-f  )f 


which  on  integration  yields 

2e^  i-e^^ 

f(£>  =  -^=57  ,  *(£)  =  ±-=37  •  (43d) 

l+e  ^  1+e  4 

It  can  be  verified  that  when  the  solution  (43d)  is  used  in  (43a)  then  the 
resulting  metric  coefficients  g^  and  g  are  equal.  Thus  the  coordinates 
£,£  are  isothermic.  The  relations  between  the  standard  spherical  polar 
coordinates  0,<J>‘and  the  coordinates  £,£  are 

0 

£  =  $  ,  £  =  in  tan  —  . 

Refer  also  to  §5.1.1. 

Example  2 :  3D  coordinates  between  a  prolate  ellipsoid  and  a  sphere. 

We  now  consider  the  case  of  coordinate  generation  between  an  inner  body 
n  =  nB  which  is  a  prolate  ellipsoid  and  an  outer  body  n  =  which  is  a 
sphere.  The  coordinates  which  vary  on  these  two  surfaces  are  £  and  £.  A 
curve  on  the  inner  surface  designated  as  £  =  £^  is 


x  =tcosh  n  cos  £  ,  y  =tsinh  n  sin  £  cos  £,  z  =rsinh  n  sin  £  sin  £. 

BO  BO  B  U 

(44a) 

Similarly  the  curve  C2  corresponding  to  £  =  £^  on  the  outer  surface  is 

nco  noo 

x  *  e  cos  £0  ,  y  ■  e  sin  £0  cos  £  ,  z  =  e  sin  £Q  sin  £  .  (44b) 


In  order  to  provide  the  solution  of  the  present  problem  with  coordinate 
contraction,  we  consider  Eqs.  (38)  and  assume 


£  *  5(£)  ,  n  *  n(n)  +  (45) 

where  £  «  0  corresponds  to  £  *  0  and  n  *  corresponds  to  n  *  nB*  Thus 
£(0)  »  0,  n ( n_)  *  0.  Under  the  transformation  (45),  the  only  nonzero 

®  Y  1  2 

components  of  are  P^  and  P22>  Writing 


M£)  -  ^  , 

<J£ 

P1  .  -  I  S! 

11  A  d£ 


6(n> 


dn 

dn 


l  de 


we  have 


(46) 
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By  taking  a  value  of  K  slightly  greater  than  one  (K  =  1.05)  we  can  have 
sufficient  contraction  in  the  n-coordinate  near  the  inner  surface.  For  the 
chosen  problem  since  the  dependence  on  £  is  simple,  we  find  that  the  generated 
coordinates  between  a  prolate  ellipsoid  and  a  sphere  are 

x  =  [AeBn^  +  C]cos  £,  y  =  DeBn^sin  £  cos  £  ,  z  =  DeBn^sin  £  sin  £ 


This  example  shows  that  the  chosen  generating  system  of  equations  (38) 
are  capable  of  providing  non-isothermic  coordinates  between  a  prolate 
ellipsoid  and  a  sphere. 


§4.  GENERATING  DIFFERENTIAL  EQUATIONS  BASED  ON  LAPLACE  EQUATIONS 

For  the  purpose  of  coordinate  generation  in  either  two  or  three  dimensions 

it  has  become  quite  popular,  particularly  after  the  publication  of  the  TTM 
2 

method  ,  to  adopt  a  system  of  inhomogeneous  Laplace'  equations  as  the 
generating  system.  The  inhomogeneous  terms  are  completely  arbitrary  and 
seemingly  there  is  no  guidance  from  the  analytical  side  as  to  how  they 
should  be  chosen.  Because  of  this  and  due  to  other  basic  reasons  it  is 
important  to  reconsider  the  formulation  of  the  problem  of  coordinate 
generation  based  on  Laplace'  system  of  equations  from  an  analytical  point 
of  view.  The  conclusions  drawn  from  these  considerations  are  that  the  set 
of  Laplace  equations 

V2x1  =  0  ,  i  =  1,  2,  3  (51) 

are  essentially  the  basis  of  the  TTM  method  rather  than  the  set  of  inhomo¬ 
geneous  equations 

vV  =  P^x1,  x2,  x3)  ,  i  =  1,  2,  3  ,  (52) 

where  P*-  are  the  specified  functions.  The  reason  for  this  conclusion 
is  that  a  coordinate  transformation  from  x*  to  any  other  system  x*,  both 
satisfying  the  same  boundary  conditions,  automatically  gives  rise  to  the  set 
of  equations  (52)  from  (51) .  Thus  as  soon  as  the  solution  of  the  system  of 
equations  (51)  under  the  constraints  of  a  body  conforming  boundary  conditions 
has  been  obtained  a  transformation 
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-i  -i.  123. 
x  =  x  (x  ,  x  ,  x  ) 

can  redistribute  these  coordinates  in  any  desired  manner. 

To  formulate  the  above  noted  ideas  analytically,  we  consider  the  formula 

12  13 

for  the  Laplacian  of  a  scalar  $  in  the  curvilinear  coordinate  system,  ' 
which  is 

v2o  =  g13(  ?  -  r*  -2J-)  .  (53) 

3x  3xJ  J  3x 

if  <j>  =  xm  is  any  curvilinear  coordinate,  then  from  (53)  we  obtain 


72*m  =  -gij  r". 


If  ♦  *  x  ,  where  x  is  any  of  the  rectangular  Cartesian  coordinate,  x,  =  x, 
mm2  1 

x.  =  y,  x,  *  z,  then  since  V  x  =  0,  we  obtain  using  (53) , 
z  3  m 

3^x  _  3x 

in  m  ,n2  r,  m 

g  J  — t — r  +  (7  x  )  — -  =  0  .  (55) 

3x  3x3  3x 

Taking  (51)  as  the  basic  generating  system,  we  get  from  (55) , 

.  .  32x 

g13  — : — =  0  .  (56) 

3x13x3 

Using  the  formulae  stated  in  §2,  we  getDx  *  0,  or 

m 


Dx  =  0  , 


Dy  =  o  , 


Dz  =  0  , 


where  the  operator  Dis  given  by 


D=  G.3..  +  G,3  +  G,3  +2G  3  +  2G-3rr  +  2G..3  „ 

l  CS  2  nn  3  4  cn  5  cc  6  nC 


(60) 
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In  two  dimensions*  g  =  1,  t—  =  0,  so  that  D  becomes 

Jj  OS 

D=  g22SCC  ‘  2gl29Cn  +  gll3nn  * 

Let  X1  be  another  coordinate  system  which  satisfies  the  same  body  con¬ 
forming  boundary  conditions  as  the  system  x1 ,  and  let 


-i  -i.  12  3  .  _  , 

x  =  x  (x  ,  x  ,  x  )  ,i  =  l,2,  3, 


(61) 


with 


det  (-^-r)  *  0 

ax3 


Then  an  analysis  similar  to  §3.2  shows  that 


3x  3x  «-£ 

_ m  _  _ m  3x 

3xj  '  as1  ax3  ' 


a2x 


m 


3  x  .-k  i.  3x  .2  l 
m  3x  3x  +  m  3  X 


3x13xj  3xk3x£  3X1  3xj  3x*  3xX3xj 


Using  the  last  expression  and  the  transformation  law 


ij  _*n  3X1  3X3 
}  =  g  ~~ 
3x  3x 


in  Eq.  (56)  we  get 


l  »  n 

_k£  m  ,  _rn  £  m 
g  “TT-T  +  9  Prn  ~J  =  0  ' 

3x 


3x  3x 


(62) 


where 


P*  - 
rn 


3xA  3x^  32x^ 


r  n  i  i 
3x  3X  3x  3x^ 


(63) 


and  is  symmetric  in  the  lower  two  indices.  If  now  in  Eq.  (55)  we  replace 
x*  by  xA,  gA^  by  gA^  and  introduce 


•Refer  to  comment  (ii)  at  the  end  of  the  paper. 


2_r  -ton 
V  x  =  -g 


mn 


=  P 


then  it  amounts  to  the  same  thing  as  taking  the  non-homogeneous  Laplace 
equations  (52)  as  the  generating  system.^  Thus  we  reach  the  conclusion  that 
essentially  Eqs.  (51)  are  the  basic  generating  equations  and  that  any  redis¬ 
tribution  of  the  solution  of  Eqs.  (51)  gives  rise  to  Eqs.  (62). 

Transferring  the  second  term  of  Eq.  (62)  to  the  right  hand  side  and  using 
the  formulae  developed  in  §2  which  are  applicable  to  all  coordinate  systems, 
we  obtain 


5xm  ■  '“A  +  G2p22  +  G3P33  +  2G4PL  +  2G5Plt3 


3x 

p  m 
2G6  P23}~ 

:IX 


(64) 


where  x  =  x,  y,  or  z,  and  n  is  the  same  operator  as  (60)  in  the  new  coordinate 
m  7,17 

system.  In  two  dimensions,  Eq.  (64)  gives  rise  to  the  familiar  forms 


Dx  -<922P11  "  2gi2P12  +  911P22)XC  “  (g22Pll  2gl2P12  +  gilP22)Xn  ' 


(65a) 


Dy  ~  ~(g22Pll  "  2gl2P12  +  gllP22)yC  *  (g22Pll  "  2gi2P12  +  gilP22)yTT  * 


(65b) 


It  must  be  noted  that  the  preceding  analysis  guides  one  to  a  proper 

i 

selection  of  the  quantities  Pr^  for  concentrating  the  coordinate  lines  in 

the  desired  regions.  This  selection,  though  still  arbitrary,  at  least 

£ 

suggests  that  the  chosen  P^  should  be  something  like  a  product  of  two  first 

and  one  second  partial  derivatives.  This  idea  is  important  in  the  adaptive 

coordinate  systems.  Furthermore,  the  preceding  analysis  also  exposes  for  the 

£ 

first  time  the  existence  of  the  cross  derivative  quantities  P  (r  f  n) 

which  do  not  appear  if  one  starts  from  the  Eqs.  (52)  and  which  may  be  important 

in  non-orthogonal  coordinates.  For  example,  in  two  dimensions  the  quantities 

P£  are 
m 


P11  ‘  +  2ECntECn  *  ‘"t1  Enn  ‘ 


PU  -  <Ei>2V  +2{rf\n  *  ,nI,2'\n  ' 


'Refer  to  comment  (iii)  at  the  end  of  the  paper. 
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P*  =  (£-)2£ff  +  2£-n-£  +  (n -)2l  , 

22  n  ££  n  n  £n  n  nn 


p5?  =  U:)2r>rr  +  2  £~  n—  Hj.  +  (n-)2^  , 

22  n  ££  n  n  £n  n  nn 


P,,  =  £-£-£r,  +  (Crn-  +  nr£-)£r  +  n-n-£  , 

12  £  n  £5  £  n  £  n  £n  £  n  nn 


p.  _  =  £r£-n,_,  +  (£-rn-  +  n^C-Jn,  +  n^n^n 
12  n  C5  (  n  £  n  £n  £  0  nn 


If  £  =  £(£)  and  n  =  n(n).  then  writing 


a  =  *£,  e  =  — 


we  get  ■ 


1  1  dA  2  1 

P11  "  '  X  -  '  P11  “  0  '  P22  "  0  ' 

a£ 


P2  =  =  I  p1  =  o  P2  =  0 

22  6  '  P12  '  12  °  ' 
dn 


which  are  exactly  the  same  as  have  been  used  in  an  earlier  paper.  In  this 


case,  writing  for  brevity 


P11  =  P  '  P2 2  ~  2 


Eqs.  (65)  simply  become 


Dx  =  -(g22PX|  +  g^Qx-)  , 


(66a) 


Dy  =  -(g22  +  91!^ 


(66b) 


1  2 

These  equations  do  not  contain  the  cross  derivative  terms  P^2>  P^2  because 
£  and  n  have  been  chosen  to  be  functions  of  £  and  n  respectively. 
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S4.1  Case  of  orthogonal  coordinates 

In  general,  for  the  generation  of  orthogonal  coordinates  it  is  not 
necessary  that  the  coordinate  functions  should  also  satisfy  the  Laplace 
equations  in  the  xyz-space.  In  this  section  after  summarizing  the  basic 
generating  equations  for  the  orthogonal  coordinates  we  have  studied  the 
effect  of  constraining  the  coordinate  functions  to  be  simultaneously  harmonic. 
The  orthogonality  conditions  are 


g .  .  =  0  for  i  ^  j  . 


(67) 


Also,  for  orthogonal  coordinates  Eqs.  (54)  simply  become 

’  ,!£l,  . 

r  dt,  h 

>3  1 


.2  1  3  hl\ 

7  n  *  7  ^  hr->  • 

va  2 


(68) 


V2,.  >1,^,  . 

JZ  H  h2 


whe  re 


' ■  SX*  +  3yy  +  ■  hl  =  ^7  •  h2  =  ^  '  h3  *  ^17  •  ^  ‘  hlh2h: 

Proceeding  directly  from  Eq.  (55)  and  using  Eqs.  (67)  and  (68)  we  obtain 


~x  =  0  ,  m  =  1 ,  2 ,  3  , 
m 


(69) 


where 


h2h3 


JL(Vi 

h. 


Note  that  the  operator 


H  and  the  Laplacian  operator 


are  related  as 


E*  =  h1h2h3  v2<i' 


for  a  scalar 
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Equations  (69)  are  those  fundamental  equations  which  every  orthogonal 
coordinate  system  must  satisfy.  A  program  of  calculation  using  Eqs.  (67)  and 
(69)  along  with  the  definitions  of  g^,  922  and  933  can  be  developed. 

§4.1.1  Case  of  orthogonal  coordinates  using  the  Laplace  equations 
Case  I :  3D  coordinates. 

If  the  generating  system  of  equations  is  taken  as 

v25  =  0  ,  v2n  =  0  ,  v2;  =  0  ,  (70) 


then  from  Eqs.  (68)  we  find  that 


hx  =  f2(C,C)f3tt,n)  ,  h2  =  f1(n,Of2(C.n)  ,  h3  =  (r,,c)  f2  (5,  5)  ,  (71) 

where  f^,  f  ,  f^  are  arbitrary  functions  of  their  arguments.  Also  the 
generating  system  (69)  for  the  Cartesian  coordinates  becomes 

32x  32  32x 

922g33  2  +  gHg33  ~2  +  911922  2  0  ,  m  ■  1,  2,  3  ,  (72) 

3^  3n  3; 

which  because  of  (71)  can  also  be  written  as 
a2  a2 

2  °X  2  °X  2  dX 

f-^0,5)  - |  +  f2<5,5)  - 1  +  f3U.n)  - f  =  0  ,  m  =  1,  2,  3  .  (73) 

35  3n  95 

Case  II:  2d  coordinates. 

For  the  case  of  2D  orthogonal  coordinates  the  equations 

v25  =  0  ,  y2p  =  0  ,  (74) 


with  the  use  of  Eqs.  (68)  yield 


922  =  a9ll  ' 


where  a  is  a  constant.  The  case  a  =  1  gives  the  corresponding  isothermic 
coordinates  which  are  conformal.  However,  by  a  straight  forward  coordinate 
transformation  of  the  isothermic  coordinates  5*0  to  another  coordinates 
5,n  we  can  have  a  coordinate  distribution  in  which  g22  f  g^.  For*  let 


5  =  5(C*n)  ,  n  =  n(5,n) 


25 


be  an  arbitrary  orthogonal  transformation.  Using  the  chain  rule  of  differ¬ 
entiation,  we  get 

C  =  CrC  +  C-n  +  ) 2  +  C--(n  )2  +  2C--C  n 

xx  c  xx  n  xx  ^  x  nn  x  Cn  x  x 

etc. ,  etc. , 

which  when  used  in  Eqs.  (74)  along  with  the  orthogonality  condition 

l  n  +  C  n  =o 
xx  y  y 


and  the  formulae 


2-  i  a  - 

VC  ^22/gll  ' 

/g  3C 


(75a) 


2-  1  3  .1-  - 

v  n  =  -T%n/g22  , 

•'g  3n 


(75b) 


<</  +  V2  ■  '  (\)2  *  <V2  =  -  5  -  ini22  • 


yield  the  equations 


■T  V2c  -  ^  <5ffe22/9n)  *  ^  <^n/522>  -  o  . 


(76a) 


^  V2n  =  _L  (n-\/g  /g  )  +  (n-Vg  /i22)  =  0  . 

35  *  ■Li  3n  n 


(76b) 


A  study  of  Eqs.  (76)  suggests  that  if  C  is  only  a  function  of  C»  and 
n  is  only  a  function  of  n,  e.g.. 


-f’ 


(C)dC  ,  n(n) 


then  Eqs.  (76a, b)  are  identically  satisfied  by  taking 


Vgu/522  =  pU)v(n) 


(76c) 


Thus 
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H 


p2(C)v2(n)522,  (76d) 


and  so  the  coordinates  £,n  are  orthogonal  but  not  conformal. 

An  important  result  from  the  preceding  analysis  is  that  if  the  orthogonal 


coordinates  are  generated  through  the  solution  of  the  Laplace  equations  (74) 


then  there  exists  an  infinity  of  transformations  £  =  £(£),  n  =  n(o)  in 


which  the  ratio  g,,/g__  is  a  product  of  a  function  of  \  and  a  function  of  n- 


This  result  is  not  in  general  true  for  coordinates  not  satisfying  the  Laplace 


equations . 


§5.  GENERATING  DIFFERENTIAL  EQUATIONS  BASED  ON  THE  RIEMANN  TENSOR 

In  any  given  space  there  are  endless  possibilities  for  the  introduction 

of  coordinate  curves.  Each  chosen  set  of  curves  determines  its  own  metric 

components.  For  example,  in  a  Cartesian  plane  besides  introducing  rectangular 

Cartesian  coordinates  x,  y,  we  also  have  endless  possibilities  for  introducing 

either  orthogonal  or  nonorthogonal  coordinate  curves.  However,  as  is  well 

known,  there  is  a  basic  differential  constraint  on  the  variations  of  g..'s 

ID 

irrespective  of  the  coordinate  system.  Since  the  curvature  of  an  Euclidean 

two-dimensional  plane  is  identically  zero,  the  basic  differential  constraint 

on  the  g. . ’s  is 
ID 

,  a  *^7  ->  *  •''57  _ 

(G-)'V,.,  =  ■—  <—  r;.>  -  £  (—  r:,)  =  o  ,  (77) 

3  1212  3n  gn  11  3£  gn  12 

where  5,n  are  any  arbitrary  coordinate  curves  in  the  plane.  Thus  no  matter 
which  coordinate  system  is  introduced  in  a  plane,  the  corresponding  matrics 

g. .  must  satisfy  Eq.  (77).  Equation  (77)  has  also  been  used  as  the  basic 

13  18 
generating  equation  for  the  generation  of  orthogonal  coordinates  in  a  plane 

12  13 

In  general,  the  Riemann  curvature  tensor  R  .  defined  as,  ' 

rjnp 

,2  a2  a2  a2„ 

,  3  g  3  g .  3  g  9  g. 

R  =  I{ _ I£E  +  _ Z2R _ llR  _  _ LIE, 

rjnp  2  3xj9xn  9xr3xP  tehy?  3xr3xn 

+  gtS(ljn,s]  [rp,t]  -  [ jp,s] [rn,t])  (78) 

defines  the  components  of  the  curvature  tensor  of  any  general  space.  If  the 
space  is  N-dimensional,  then  the  number  of  components  R  .  are  given  by 


Thus  for  N  *  2  there  is  one  distinct  surviving  component  stated  in  Eq.  (77). 
However,  for  N  =  3,  it  has  six  distinct  components 


**1212 '  **1313 '  **2323'  **1213'  **1232'  **1323 


If  the  3D-space  is  Euclidean,  then  its  curvature  is  zero,  so  that  the  six 
equations 

**1212  =  °  '  **1313  =  °  '  **2323  =  °  ' 

**1213  =  °  '  **1232  =  °  '  **1323  =  ° 


(79) 


determine  the  differential  constraints  for  the  six  metric  coefficients  g..  in 

ID 

any  coordinate  system  introduced  in  an  Euclidean  space.  These  equations  in 


**1212 


1313 


2323 


led  form 

are  as 

follows : 

.2 

.2 

.2 

3  gll 

3  gl2 

3  g22 

3n2 

Z 

3?3n 

+  H2 

.2 

.2 

.2 

3  gll 

3  g13 

3  g33 

3C2 

z 

3C3C 

V 

_2 

.2 

.2 

3  g22 

3  g23 

3  g33 

+  2g  S  ( [ 22 , s ] [ 11 , t ]  -  [ 12  ,  s ] [ 12 , t ] )  =  0  , 


(80a 


+  ([33, s] [ 11 , t ]  -  [ 13 ,s ] [ 13 , t ] )  =  0  , 


(80b 


3C 


3n3; 


+  2g  S  ( [ 33 , s ] [ 22 , t]  -  [ 23 ,s ] [ 23 ,t ] )  =  0  , 


3n 


(80c 


**12: 


32gH 

32gl2 

a2g 

g13 

32g23 

3n3; 

3C3  C 

3£3n  + 

3C2 

32g22 

32g12 

32g23 

32g13 

3C3C 

3n3s 

3£3n  + 

3n2 

92g 

g33 

a2g 

g13 

a2g 

g23 

— — - —  + 

(80d 


+  2g  S ( [ 22 ,s ] [ 13 , t ]  -  [23, s] [12, t]) 


1323  3C3n  3n3C  3£3;  -r2 


(80e 


+  2g  S([33,s)[12,t]  -  [23 ,s ] [ 13 ,t] )  = 


(80f 


where  [ij,k]  are  the  Christoffel  symbols  of  the  first  kind  defined  in  (8a). 

Equations  (80)  are  those  consistent  set  of  partial  differential  equations 
which  must  always  be  satisfied  by  the  metric  coefficients  g. ..  In  the  3D  case 
Eqs.  (80)  are  six  equations  in  six  unknowns  and,  therefore,  they  form  a  closed 
system  of  equations.  In  contrast,  for  the  2D  case  there  is  only  one  equation 
(Eq.  (77))  and  three  unknowns  g11f  g12,  g22  and  therefore  some  constraints 
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are  needed  to  turn  Eq.  (77)  (such  as  orthogonality ^8)  into  a  solvable 
equation.  This  author  is  not  aware  of  any  numerical  solution  of  the  complete 
set  of  equations  (80) ,  though  there  are  some  possibilities  of  developing 
solution  algorithms  using  Eqs.  (80)  as  the  core  equations.  For  example, 
in  the  problem  of  obtaining  the  3D  coordinates  for  the  configuration  of 
Fig.  1,  one  can  judiciously  choose  g^,  ^13'  ^33  kased  on  t^e  given  bound¬ 

ary  data  for  the  whole  field  and  then  solve  Eqs.  (80)  for  the  remaining 
coefficients  g^/  ^23'  ^12“  It:  s^ou^d  also  be  noted  that  in  any  physical 

problem,  e.g. ,  the  Navier-Stokes  problem,  one  only  needs  the  metric  coeffi¬ 
cients  and  their  derivatives  (Christoffel  symbols)/  which  become  available 
after  solving  Eqs.  (80).  Nevertheless,  for  graphical  and  other  purposes, 
one  also  needs  the  functions  x(£,n,0  etc. 

To  obtain  the  Cartesian  coordinates  on  the  basis  of  the  available  g. .  *s, 

ij 

we  introduce  the  unit  base  vectors  A .  as 


A.  =  a./vg. .  ,  no  sum  on  i. 

li 


Let  the  components  of  A^  along  the  rectangular  Cartesian  axes  be  denoted  as 

u. ,  v. ,  w. ,  so  that 
X  X  X 


where 


*i  =  (V  V  V  ' 


ui =  v/gn '  vi =  v/gu '  wi =  v>gn  ' 


uo  =  »  wo  =  zJ 


n'  y22  ' 


U3  =  X^//g33  '  V3  =  yC//g22  '  W3  =  * 


Knowing  u^,  v^,  w^,  it  is  possible  to  evaluate  the  Cartesian  coordinates 
through  the  line  integrals 


r  =  (A^i 


911  d*  +  *2^22  dn  +  Vg33  d;)  * 


The  determination  of  u^,  v^,  w^  is  a  separate  problem  which  we  now 
consider.  First  of  all  using  (81)  in  Eq.  (8c),  we  get  a  system  of  first 
order  partial  differential  equations 


where,  as  before,  there  is  no  sum  on  the  repeated  index  i.  Equations  (84) 

form  a  system  of  27  first  order  PDE's  in  nine  independent  variables  u^, 

Vi*  Wi'  T^s  system  of  equations  is  overdetermined  and  thus  its  solvability 

should  depend  on  certain  compatibility  conditions.  According  to  a  theorem 

19 

on  the  overdetermined  system  of  equations  ,  if  the  compatibility  conditions 
hold  then  the  solution  of  Eqs.  (84)  exists  and  is  unique.  The  conditions 


3xm3x:) 


3x3  3xm 


for  all  values  of  i,  m,  and  j  are  the  compatibility  conditions.  To  prove  (85) 
we  use  Eq.  (8c) ,  which  on  cross  differentiation  yields 

i2  i2 

3a.  3a.  { 

- - T*i-  -  K  .  a.,  (86) 

axV  ax’a*’  •1")! 

l  12 

where  R  .  .  is  the  Riemann-Christof fel  curvature  tensor  and  is  related 

.imj  . 


with  the  Riemann's  tensor  R. 

ijki 


Evidently  in  our  present  case  R  .  .  =  0, 

.  lm  j 


since  the  space  is  Euclidean.  Inserting  (81)  in  (86)  we  find  that  Eq.  (85) 
are  identically  satisfied. 

It  is  interesting  to  note  that  for  a  two-dimensional  curvilinear  coordinate 
system  there  is  no  need  to  solve  the  system  of  equations  such  as  (84) .  in 
this  case  the  single  differential  equation  with  =  g 

a  6  r?,  ,  ■'«  rj, 

■ma  ■  -  0 

implies  the  existence  of  a  single  function  a(£,n)  such  that 


a£  ~  ?! 


Consequently 


u^  ■  cos  a,  v^  =  -sin  a,  u^  ■  cos(a-0)  ,  v2  ■  -sin(o-0)  , 


3 


where  a  is  the  angle  made  by  the  tanget  to  the  coordinate  line  n  *  const,  in 
a  clockwise  sense  with  the  x-axis ,  and 


cos  e  =  g12/^ng22 


is  known. 


§5.1  Case  of  orthogonal  coordinates 


We  again  return  to  the  case  of  3D  orthogonal  coordinates.  Refer  also  to 
§4.1.  Under  the  constraint  of  orthogonality. 


gl2  *  913  =  g23  =  0  '  [12'3)  =  (13'2]  =  [23'1]  =  0  ' 


ri2  P13  *  r23  ~  °  '  9  =  9lig22g33  ’ 


the  set  of  equations  (80)  reduce  somewhat.  They  are 


33/g,.g 


1133 


22  ^ll^ 


11.  22  33 


2  on  g 


3C  30 


35on  2  g^  an 

which  are  the  Lame's  equations 
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§5.1.1  The  case  of  isothermic  coordinates. 

Isothermic  coordinates*  in  a  surface  embedded  in  a  3D  Euclidean  space 
are  those  coordinates  in  which  the  metric  coefficients  g^  and  g^3  in  the  sur¬ 
face  n  =  const,  are  equal.  That  is,  the  element  of  length  ds  on  n  =  const.  : 
given  by 


(ds,2„==onst.  '  «uKd{)2  *  ■ 


where  are  chosen  to  be  the  surface  coordinates.  Setting 


g33  =  gn,  and  g22  =  F(n) 


in  Eqs.  (88),  we  obtain  the  basic  equations  for  g  ,  which  are 


a  ,  i  3gn.  .  3,1  3gn.  i  ,3gnx2  =  0  , 

— (- - —)  +  —  (“ - — )  +  — — (-T— ) 


gn  H  3;  9U  a;  2pgn'  3r) 


(89 


li 


_i(_J - 

3n  3n 


)  =  0  , 


11 


-  0  . 


3C  gn  3n 


3  ,  1  9gil 


H  gn  3n 


)  -  o  . 


It  can  easily  be  verified  that  the  only  solution  of  Eqs.  (89c, d)  is 


gi;L  =  [a+P(n>]  f(5*C)  ,  a  *  const. 


Thus  from  (89b) 


P(n)  -  (~)  ■ 

dn 


Substituting  (90)  and  (91)  in  Eq.  (89a)  ,  the  differential  equation  for 


(89 


(89 


(89 


(90 


(91 


*Refer  to  the  comment  (iv)  at  the  end  of  the  paper. 
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f(5,5)  becomes 


J_(I  ii)  +  _L(I  M} 

35  f  35  35  f  35 


+  2f  =  0 


(92) 


14 

In  Kreyszig  we  have  the  result  that  if  in  a  portion  of  a  surface 
isothermic  coordinates  can  be  introduced  then  that  portion  of  the  surface  can 
conformally  be  mapped  onto  a  plane.  Thus  in  effect  the  solution  of  Eq.  (92) 
provides  that  mapping  function  which  conformally  maps  a  surface  onto  a  plane. 
As  a  verification  of  the  above  conclusion,  we  verify  that  the  function 


4e 

(1+e 


25 


(93) 


is  a  solution  of  Eq.  (92) .  This  function  is  related  with  the  isothermic 
coordinates  on  a  sphere.  Using  the  parametric  equation  of  a  sphere 


x  =  [a+P(n)J  cos  0,  y  =  (a+P(r,)]sin  9  sin  z  =  [a+P(n)]sin  6  cos  4> 


and  writing 


5  =  ?  ,  5  =  £n  tan  -  , 


where  0  <  <p  <  2 71  and  0  <  6  <  it,  we  obtain 


g33  ~  911 


4  (a+P)2e2^ 

(l+e2S2 


Thus  the  equations 


x  = 


(a+P) (l-e2;) 


1+e 


25 


y  = 


2  (a+P)e^sin  5 


1+e 


25 


z  = 


2 (a+P) e^cos  5 


1+e 


25 


(94) 


represent  a  sphere  of  radius  a+p(n)  in  terms  of  the  isothermic  coordinates 
5 #5  in  the  surface.  Since  P ( n)  is  an  arbitrary  function  of  r\»  we  have  the 
capability  of  prescribing  a  suitable  function  P(n)  to  achieve  any  sort  of 
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contraction  or  expansion  in  the  field.  It  looks  that  the  representation  (94) 
should  prove  useful  in  the  computational  problems  associated  with  a  sphere. 

Comment  (i) ; 

As  a  further  justification  for  the  consistency  of  the  set  of  Eqs.  (29)  it 
has  been  shown  below  that  these  equations  can  be  combined  to  obtain  the  equation 
for  a  surface  z  =  z(x,y)  in  the  well  known  form 

az  -  2gz  +  y z  =  2HM  ,  (i) 

xx  xy  yy 

where 

2  2 

2H  =  =  R/G3  ,  M  =  1+p  +q  ,  p  =  z^,  q  =  z^  , 

a  ■  (l+q2)/v1w  ,  g  =  pq/^M  ,  y  =  (l+p2)/*'?!  . 


First  note  the  following  definitions  and  identities : 

G3  =  9lig22  '  (gi2}  '  X  =  '  Y  =  •  Z  =  ' 


(x,x)  =  (1-X2)G3  ,  A1<x,y)  =  -XYG3  ,  A  (y,y)  =  (1-Y2)G3  , 


where 


Va'b)  =  g22a*  -  gi2(acbn  +  anV  +  gnanbn  *  JJ 

Calculating  z^,  z^n'  znn  from  z^*  ZT]>  substituting  these  expressions  in 
Eq.  (29c)  while  using  the  equations  in  (ii)  and  Eqs.  (29a, b)  we  recover 
Eq.  (i)  given  above. 

We  now  compare  the  equations  obtained  by  Thomas6  with  those  of  Warsi11. 
Thomas'  equations  in  the  present  notation  are 

£x  +  2pG  H/i^T  =  0  ,  £y  +  2qG.H/<^M  =  0  ,  where  G,  =  (x_y  -x  yr)2M  ,  (iii) 
3  3  3  (  11  H 


which  are  exactly  the  same  as  Eqs.  (29a, b)  of  this  paper.  It  must,  however, 

i 

be  pointed  out  that  the  derivation  of  Eqs.  (iii)  involves  fo>’r  steps: 

(a)  orthogonality  of  (  with  £,n,  (b)  vanishing  of  the  principal  curvature  of 
(-lines,  (c)  elimination  of  an  arbitrary  parameter  (which  may  be  zero), 

(d)  prescription  of  z(x,y)  for  mh&  surface  to  be  generated. 
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Comment  (ii) ; 

In  two  dimensions  another  differential  system  is  provided  by  a  first  order 
20 

Beltrami  equation  ,  which  in  the  complex  form  is 


f-  -  H(z,z)fz  =  0  , 

(i) 

where 

f  =  f(z,z)  , 

z  =  x+iy,z  =  x-iy,i  =  S-l  . 

Writing 

f(z,z)  =  £ (x,y )  +  in(x,y)  ;  H(z,z)  =  p(x,y)  +  iv(x,y)  , 

(ii) 

we  obtain 

the  following  two  real  equations  from  (i) : 

'  8-x  *  Y5y  • 

(iii) 

\  '  “Cx  *  65y  ' 

iiv) 

where 


a  =  [  (1—  p) 2  +  v2]/ A  ,  8  *  -2v/A  ,  y  =  [(1+u)2  +  v2]/A  ,  A  =  l-(y2  +  v2)  . 
Note  that 

ay  -  B2  =  1  , 
a  +  y  =  2 (2-AJ/A 

A  quasiconformal  mapping  becomes  conformal  when  H  =  0,  or  equivalently 
a  =  y  =  1#  B  =  0.  The  resulting  equations  are  then  the  Cauchy- Riemann 
equations 


and  then  f<z)  is  an  analytic  function  in  the  domain  D. 

Equations  (iii)  and  (iv)  can  be  inverted  so  that  only  the  partial 
derivatives  of  x  and  y  appear.  Thus 


4.  • 
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x  =  Bx  -  ay  ,  (v) 

n  £  £ 

y  =  yx  -  By.  •  (vi) 

n  £  £ 

For  Eqs .  (v)  and  (vi)  it  is  important  to  write  a,  g,  y  in  terms  of  the  metric 
coefficients,  which  are^'^2, 


A  =  4  /g/  [  2  /g  +  (gn  +  g^)]  , 

a  +  y  =  (gn  +  g22)/*/g  , 

<**Y  =  [glx  +  g22  +  t  <9X1  +  922)2  -  4(1  +  B2)g}>5]/2»/g"  . 


Comment  (iij) ; 

As  is  expected,  Eq.  (82)  can  be  reduced  to  the  form 


-k£ 

g 


32x 


m 


,-k._£ 
dx  3x 


(72lr) 


3x 

m 


0 


by  using  the  formula 


-I  „-q  _r 

3  X  _  J.p  3x  rp£  bX  3x 

sx^x5  =  ij  Sx*  "  qr  3X1  3xj 


in  the  expression  for 


32x 


m 


However,  for  gaining  a  new  insight  into  the 


structure  of  the  redistribution  terms  it  looks  profitable  to  keep  the  form 
(62)  with  defined  in  (63) . 


Comment  (iv) : 

Generation  of  isothermic  coordinates  can  also  be  achieved  by  the  method 

~1  '2 

detailed  in  Ref.  14.  Let  x  and  x  be  some  sort  of  coordinates  introduced 
in  a  portion  of  the  surface  (for  example  from  the  subroutine  developed  by 

2i  12 

Craidon  ),  and  let  x  ,  x  be  the  desired  isothermic  coordinates.  Then 


i 

x 


i/'l  ~2s 

X  (x  ,x  ) 


Because  of 


being  isothermic, 


we  have 


g22  =  911  * 


36 


Using  the  transformation  law  for  the  covariant  and  the  contravariant  metric 
tensor  components  we  get 


(i) 


where 


'11 


=  0 


'22 


=  0 


'12 


=  +1 


'21 


=  -1 


(ii) 


From  Eqs.  (i)  and  (ii)  we  find  the  second  order  differential  equations 


3  ,  rr  ~i j  3x\  .  ... 

— 7-  (/g  g  — T-)  =  0  ,  (xx 

3x  3*3 

where  k  =  1,  2.  Note  that  in  the  Eqs.  (i)  -  (iii)  the  indices  range  over  the 
values  1,  2. 

Equations  (iii)  provide  two  linear  uncoupled  equations  for  the  deter¬ 
mination  of  the  isothermic  coordinates,  since  the  values  g1"1  of  g1"*  are 
known  a’  priori. 


CONCLUSIONS 

Three  distinct  methods  of  numerical  coordinate  generation  based  on  PDE's 
have  been  analyzed  in  detail .  In  the  two  newly  proposed  methods ,  viz . ,  the 
methods  discussed  in  S § 3  and  5,  some  useful  results  have  been  obtained  by 
looking  at  the  generating  system  of  equations  as  a  system  of  forcing  differ¬ 
ential  relations  among  the  metric  coefficients  g^ .  For  example,  in  the 
method  of  S3  and  g_'s  are  forced  to  satisfy  Eqs.  (27)  (refer  also  to  their 
forms  in  Eqs.  (20)).  In  the  method  of  §5,  the  g_'s  naturally  satisfy  Eqs.  (80) 
since  the  space  is  intrinsically  Euclidean.  In  the  TTM  method  discussed  in 
§4  the  generating  Laplace  or  Poisson  equations  also  amount  to  specifying  a  set 
of  differential  constraints  on  the  g^'s. 

In  the  process  of  obtaining  the  above  noted  results  a  number  of  other 
results  and  equations  have  been  obtained  which  should  be  satisfied  by  all 


coordinate  systems.  For  example,  the  orthogonal  coordinates  in  an  Euclidean 
space  must  satisfy  Egs.  (69),  (88),  and  the  nonorthogonal  coordinates  must 
satisfy  Eqs.  (80) ,  no  matter  which  method  is  used  to  generate  them.  In 
effect  all  these  results  provide  enough  material  for  proposing  more  efficient 
calculation  algorithms  for  the  coordinate  generation  on  a  computer. 
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NUMERICAL  GENERATION  OF  THREE-DIMENSIONAL  COORDINATES  BETWEEN  BODIES  OF 
ARBITRARY  SHAPES 

Z.  U.  A.  WARSI*  AND  J.  P.  ZIEBARTH^ 

Department  of  Aerospace  Engineering,  Mississippi  State  University, 

Mississippi  State,  Mississippi  39762,  USA 

INTRODUCTION 

This  paper  is  devoted  to  the  numerical  solution  of  a  set  of  second  order 
elliptic  partial  differential  equations  for  the  generation  of  three-dimensional 
curvilinear  coordinates  between  two  arbitrary  shaped  bodies.  The  central 
idea  of  the  method  is  to  generate  a  series  of  surfaces  between  the  given 
inner  and  the  outer  boundary  surfaces  and  then  to  connect  these  surfaces  in 
such  a  manner  so  as  to  have  a  sufficiently  differentiable  three-dimensional 
coordinate  net  in  the  enclosed  region. 

The  basic  analytical  foundation  of  the  present  method  has  already  been 
laid  out  by  Warsi  in  52  of  Ref.  1.  However,  it  is  important  to  state  here 
that  the  proposed  equations  for  the  numerical  solution  form  a  consistent 
set  of  second  order  elliptic  equations  which  are  a  consequence  of  the 
equations  of  Gauss2  for  a  surface.  Additional  constraints  are  then  imposed 
which,  besides  yielding  the  simplest  form  of  equations  for  numerical 
purposes,  also  preserve  the  essential  geometric  properties  of  the  generated 
surfaces. 

Formulation  of  the  mathematical  model 

To  fix  ideas,  let  it  be  desired  to  generate  the  coordinates  between  the 
two  surfaces  designated  as  n  "  Hg  (the  inner  surface)  and  n  ■  (the  outer 
surface)  respectively  as  is  shown  in  Fig.  1  .  The  two  coordinates  which 
vary  in  these  two  surfaces  are  then  labeled  as  £  or  I  and  £  or  K.  The 
surfaces  n  *=  nn  and  n  ■  n  are  the  known  surfaces  in  which  the  Cartesian 
coordinates  r  *  (x,y,z)  are  given  as  functions  of  £  and  £,  that  is, 

r  *  r0(£,£)  *  r  B  r00(C,C) 

are  known  either  numerically  or  analytically.  The  method  to  be  discussed 
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generates  a  surface  for  each  fixed  value  of  C  or  K  starting  from  a  curve  on 
B  and  ending  at  the  corresponding  value  of  (  or  K  on  the  outer  boundary 
surface.  Refer  to  Fig.  lb. 


Figure  L(a)  .  Coordinates  £  and  <;  on  the  given  surfaces,  (b)  the  generated 
surface  £  =  const. 

Referring  to  Eq.  (18)  in  Warsi1,  we  now  impose  the  restrictions 


A25  =  G3<2912T12  "  922T11  "  911T22)  =  °  ' 


(1) 


A2n  "  G3{2912T12  "  922T11  _  911T22)  *  °  ' 


(2) 


for  t  =  const.  In  Eqs.  (1)  and  (2)  A  is  the  Beltrami's  second  order 

12  2 
differential  operator  '  , 


and 


3 


T2  -  —  (a  -^i+a  f"22  2  ffl£n 

T22  2G  l9ll  3n  912(  H  "  2  Srj  }1, 


T1  _  J_  .  !fl2  !f22,  *22. 

T22  2G3  922  2  3n  "  3C  ”  912  9n  1  ' 


T2  -_L_fa  o!!ii  Ohh  *u. 

T11  2G  lgll(2  SC  Sn  912  SC  ' 


.1  .  _s_(a  !!ii  ^22, 

12  2G3  922  Sn  "  912  SC  '  ' 


r2  =  — i — (g  Zli 

12  2G,  yll  sc 


’12  Sn 


Based  on  the  structure  of  the  Christoffel  symbols  in  Eqs.  (4)  we  conclude 

PY 

that  the  constraining  equations  (1)  and  (2)  are  essentially  a  set  of  differen¬ 
tial  constraints  on  the  variations  of  the  metric  coefficients  g  .  Thus 

Gp 

under  the  constraining  equations  (1)  and  (2),  the  three  equations  for  the 
generation  of  the  Cartesian  coordinates  x,  yf  z  can  be  obtained.  Below 
we  write  the  equations  when  it  is  desired  to  have  a  concentration  or  expansion 
in  the  coordinates  C  and  n,  (refer  to  Eqs.  (38)  in  Warsi1) .  For  brevity  of 
notation  we  use  the  same  coordinates  (C*n>C)  either  with  or  without 
coordinate  redistributions.  The  equations  are 

£x  =  XR  ,  (5a) 


£y  =  yr  , 


where 


£z  =  ZR  , 


£-  «22S«  -  +  ' 
*  *  Vn  •  V«)//57  ' 

Y  -  (x^  -  x^l/YSJ  , 


Z  -  (x5yn  -  Xny?)/-S7 


9 


(7c) 


4 


*  -  «*4  *  »yc  -  22;)  (gnr^  -  2Sl2rJ2  +  g22r^)  , 


r3  =  _L,G  !!ii  .  r  (2  Uhl  agll.  .  r  (  2  !!l3  *11.. 

rn  2g  G5  g6(2  a?  ■  an  G3(  2  as  "  as  )J  '  (9a) 


r3  -  _i.fr  9911  .  r  9g22  3g13  9g23  5gl2,  , 

ri2  "  2g  G5  3n  G6  as  G3  3n  +  3S  “  3C  ' 


3  1.  .  9g!2  3g22  9g22  3g23  3g22, , 

r22  "  2g  G5  2  3n  ‘  3S“^  +  G6  3n  +  G3(2  3n  "  3?  )J  ' 


gl2923  “  913g22  ' 


(10a) 


G6  "  912gl3  “  911923  ' 


(10b) 


9  =  933G3  +  913G5  +  g23G6  ' 


(10c) 


and  G3  has  already  been  defined  earlier  in  (3d) . 

A  successful  program  of  calculations  based  on  the  set  of  Eqs.  (5)  -  (10) 
now  rests  on  how  effectively  one  can  devise  a  calculation  method  for  the  first 
partial  derivatives  r^  =  (x^,  y^,  z in  the  field.  In  this  connection  we 

first  note  that  based  on  the  prescribed  values  r  (S,C)»  r  (£»G)  the  partial 

,, 00 

derivatives  with  respect  to  S  and  £  of  any  order  can  be  evaluated  on  the 
given  bodies.  Thus  we  must  somehow  connect  the  evaluation  of  r^  in  the  field 
with  the  partial  derivatives  in  the  surface.  To  maintain  the  intrinsic 
geometrical  properties  of  the  C-lines  in  the  field  with  the  c-lines  of  the 
inner  and  the  outer  boundaries,  we  consider  the  differential  equations  for 
the  surfaces  n  *  const.  Following  the  method  in  Warsi*'3  we  find  that  the 
coordinates  £,  t,  in  any  surface  (including  the  given  boundaries)  must  satisfy 
the  equations 

933*CC  ”  29135u  +  911?«  +  <G2A2  “  G2(kl  +k2  '  (11) 

where  the  enclosed  superscript  (2)  in  Eq.  (11)  means  that  all  the  quantities 
have  been  evaluated  on  the  surface  n  ■  const.  Also 


G2  “  9U933  '  (g13)  * 


(12a) 
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g2 (k1+k2)  -  g33u(2)  -  2g13T(2>  +  gns(2)  ,  (12b) 


U(2)  -n<2>.ru,T(2>=n(2,-ru,S(2)=n(2).ru 


n(2)  ,  /y<2>  v(2)  ?(2)1 

n  *  (X  ,  Y  ,  Z  )  , 


where 


X  *  (y;zC  ‘  yCZC>/*^2  ' 


(14a) 


Y  2)  =  <x^  -  x^)/*^  » 


(14b) 


z  2  =  (x;yc  -  x?y?)/»^  . 


(14c) 


It  may  be  noted  that  (kj+k2)/2  is  the  mean  curvature  and  S,  T,  U  are  the 
coefficients  of  the  second  fundamental  form  of  the  surface  n  =  const. 

Based  on  Eq.  (11),  we  now  formulate  the  following  weighted  integral 
formula  for  the  evaluation  of  r^  in  the  field. 


!c  =  +  f2(n>(!«)“]  dC  ' 


(15a) 


where 


<r„>.  -  1-^  <k‘2V2>>n,2> 


B,®  lgxl  1 


911  gll  “ ^ 


42)U*c]b-' 

gn  2  _C  B,» 


(15b) 


fl<nB>  =  1  '  fl<n»)  *  0  *  f2(nB>  *  °  '  f2<n®)  =  1  * 


(15c) 


The  functions  f^(n)  and  f ^ ( n)  must  satisfy  the  conditions  (15c)  and  should 

be  chosen  to  reflect  the  effect  of  the  coordinate  redistribution  function  Q 

appearing  in  Eq.  (6) .  It  is  also  to  be  noted  that  the  coordinate  C  need  not 

(2) 

in  general  satisfy  the  Beltrami  equation.  That  is,  in  general  CKO. 
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Numerical  solution  of  the  equations 

The  numerical  method  used  in  this  research  for  solving  the  system  of 
Eqs.  (5)  is  the  method  of  finite  difference  using  the  point-SOR.  First  the 
coordinates  5  and  4  for  both  the  inner  surface  (n  =  r^)  and  the  outer  surface 
tn  *  I*)  are  to  be  generated  using  the  available  x,  y,  z  values  for  these 
surfaces  either  analytically  or  by  a  computer  program  developed  by  Craidon.4 
In  this  research  we  have  used  both  the  analytical  methods  where  possible, 
and  also  the  subroutine  in  Ref.  4  to  generate  the  given  body  surface  coordin¬ 
ates,  with  equal  success.  Three  practical  problems  have  to  be  resolved  before 
an  effective  solution  algorithm  for  Eqs.  (5)  can  be  developed.  They  ares 

(i)  a  specification  of  the  functions  f^(n)  and  f2(n)  appearing  in  Eqs.  (15), 

(ii)  specification  of  the  redistribution  functions  (concentration  or  expansion 
functions)  P  and  Q,  and  (iii)  a  method  to  obtain  the  same  coordinates  on  the 
inner  and  outer  boundaries.  We  now  discuss  each  problem  in  succession. 

(i)  Before  discussing  the  specifications  of  f^(n)  and  f^ln)  we  may  state 
that  each  value  like  n  *  n0  and  n  =  nto  is  a  parameter  to  start  with  rather 
than  an  integer.  The  difference  n,,,  ~  nB  is  the  most  important  difference 
and  is  known  as  the  "modulus  of  the  domain."  The  determination  of  n  “  is 

00  3 

a  formidable  problem  in  three  dimensions  but  fortunately  there  is  no  need 
for  it  in  the  case  of  numerical  coordinate  generation.  Writing 

n -n 

Z  =  -  ,  (16a) 

we  find  that  the  function  f^  defined  in  (15c)  should  be  a  function  of  Z  only, 
so  that 


fx  (1)  =  1  ,  f  (0)  -  0  , 


(16b) 


and 


f2(Z)  -  1  -  fx(Z)  .  (16c) 

In  the  present  computations  we  have  taken  f^  and  f 2  as  linear  functiorsof  Z, 
that  is 

f2(Z)  =  Z.  (17a) 


Other  simple  possibilities  which  have  been  tried  are 


Though  for  convex  surfaces  the  method  of  spherical  projection  seems  to  be 

most  desirable,  we  have  for  the  present  investigations,  used  the  geometrical 

method  of  first  surrounding  the  inner  body  by  a  sphere  of  diameter  equal 

to  the  major  length  of  the  inner  body.  Next,  each  point  (x,y  ,  z  )  on  the 

B  B 

inner  body  is  projected  to  a  point  (x,y  ,  z  )  on  the  sphere  surrounding  the 

s  s 

inner  body.  The  correspondence  between  the  inner  and  outer  body  is 

then  made  by  extending  a  straight  line  from  the  center  through  (x,y  ,  z  ) 

s  s 

to  a  point  (x  ,  y  ,  z  )  on  the  outer  sphere. 

A  number  of  program  runs  have  been  made  for  prolate  ellipsoids  of  various 
thicknesses  surrounded  by  sphere  of  large  radii.  Also  a  thin  body  of 
revolution  with  circular  sections,  resembling  the  fuselage  of  an  airplane, 
surrounded  by  a  sphere  has  been  considered.  These  numerical  results  with 
and  without  coordinate  concentration  are  shown  in  Figs.  2-7. 


*■11 


-1  0  -tl 


-1  0  +1 


Figure  2.  Inner  body  a  thick  prolate  ellipsoid  with  major  axis  2  and 
minor  axis  /J  surrounded  by  a  sphere  of  radius  4.  (a)  Coordinate  contours 

for  a  section  £  *  const.  (K  =  11)  for  all  (C,n)  or  (l,J)  values, 

(b)  for  a  section  £  =  const.  (I  =  1)  for  all  (n,C)  or  (J,K)  values.  In 
both  cases  no  contraction  in  n»  *  *  1. 


Figure  3.  Cases  (a)  and  (b)  of  Fig.  2,  with  contraction  in  n»  K  =  1.05. 
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Figure  4.  Inner  body  a  thin  prolate  ellipsoid  with  major  axis  1.02, 
minor  axis  0.201  surrounded  by  a  sphere  of  radius  1.5.  (a)  Coordinate 

contours  for  a  section  c  *  const.  (K  *  11)  for  all  (£,n)  or  (I,J)  values 
(b)  for  a  section  £  =  const.  (I  *  1)  for  all  (n,0  or  (J,K)  values.  In 
both  cases  no  contraction  in  n,  K  =  1. 


K  =  1.02. 
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Figure  5.  Cases  (a)  and  (b)  of  Fig.  4,  with  contraction  in  n» 


Figure  6.  Inner  body  a  thin  body  of  revolution  with  circular  sections 
having  major  axis  2  and  minor  axis  0.1313  surrounded  by  a  sphere  of  radius 
2.5.  (a)  Coordinate  contours  for  a  section  5  *  const.  (K  =  15)  for  all 

(5,n)  or  (I,J)  values,  (b)  for  a  section  f.  *  const.  (I  *  1)  for  all  (n»t) 
or  (J,K)  values.  In  both  cases  no  contraction  in  n»  K  *  1. 
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Figure  7.  Cases  (a)  and  (b)  of  Fig.  6,  with  contraction  in  n#  <  =  1.005. 

CONCLUSIONS 

This  paper  has  been  devoted  to  the  numerical  solution  of  a  set  of  elliptic 
equations  for  the  purpose  of  numerically  evolving  a  series  of  surfaces 
and  the  intersecting  surfaces  in  arbitrary  three-dimensional  regions  in 
space.  The  most  difficult  part  of  such  a  program  is  the  generation  of 
surfaces  between  any  two  given  surfaces.  This  has  been  considered  here 
for  thick  and  thin  prolate  ellipsoids  and  a  body  of  revolution  forming 
the  inner  bodies  and  a  sphere  forming  the  outer  boundary.  Many  successful 
numerical  algorithms  can  be  developed  using  the  proposed  equations  as 
the  core  equations  for  providing  the  coordinates  around  a  complete  aircraft 
and  other  aerodynamical  shapes. 
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